package fem2.material.cohesive_law2;

import fem2.Element;
import fem2.Variable;
import fem2.element.crack_sda2.SDAElement;
import fem2.enu.VariableName;
import fem2.jlinalg.DiagonalMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.MathUtilities;
import math2.MutableDouble;

/**
 * Mode I cohesive law based on Sancho paper: An embedded crack model for finite
 * element analysis of concrete fracture.<br>
 * Note: this class works with constant opening only
 * 
 * @author hbui
 * 
 */
public abstract class SanchoCohesiveLawImplex extends CohesiveLaw {

	public SanchoCohesiveLawImplex(double Gf, double ft) {
		super(Gf, ft);
	}

	/**
	 * compute value of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLaw(double x);

	/**
	 * compute value of derivative of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLawDerivative(double x);

	@Override
	public IVector getJumpTraction(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		IVector w = getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP);
		IVector t = (IVector) w.clone();

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);
		double wt_alg = getVariable(e).get(k)
				.getScalarVariable(VariableName.EQUIVALENT_JUMP_ALGORITHMIC).getValue();
		if (wtilde < MathUtilities.tol6 || wt_alg < MathUtilities.tol6) {
			double ft = computeSofteningLaw(0);
			t = new ArrayVector(dim);
			t.set(0, ft);
			return t;
		}

		double coeff = computeSofteningLaw(wt_alg) / wt_alg;
		BLAM.scale(coeff, t);

		return t;
	}

	@Override
	public IMatrix getJumpStiffness(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		IVector w = getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP);
		IMatrix C;

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);
		double wt_alg = getVariable(e).get(k)
				.getScalarVariable(VariableName.EQUIVALENT_JUMP_ALGORITHMIC).getValue();
		if (wtilde < MathUtilities.tol6 || wt_alg < MathUtilities.tol6) {
			C = new DiagonalMatrix(dim);
			BLAM.zero(C);
			return C;
		}

		C = new DiagonalMatrix(dim);
		double coeff = computeSofteningLaw(wt_alg) / wt_alg;
		DiagonalMatrix.set(C, coeff);

		// Debugger.watch(e.getId() + ", wt_old (t) = " + wt_old);

		return C;
	}

	@Override
	public void Initialize(Element e) {
		int ngc = ((SDAElement) e).getNumGaussPointOnCrack();
		int dim = e.getMeshPart().getDimension();

		ArrayList<Variable> vars = new ArrayList<Variable>();

		for (int i = 0; i < ngc; i++) {
			Variable v = new Variable(2);

			v.registerVectorVariable(VariableName.DISCONTINUITY_JUMP, dim, false);
			v.registerScalarVariable(VariableName.EQUIVALENT_JUMP, true);
			v.registerScalarVariable(VariableName.EQUIVALENT_JUMP_ALGORITHMIC, false);

			vars.add(v);
		}

		setVariable(e, vars);
	}

	@Override
	public void InitializeSolutionStep(Element e) {
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		int ngc = ((SDAElement) e).getNumGaussPointOnCrack();
		IVector w = ((SDAElement) e).getCrackOpeningParameter();
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);

		for (int k = 0; k < ngc; k++) {

			/*
			 * save the discontinuity jump
			 */
			BLAM.copy(w, getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP));

			/*
			 * explicit update
			 */
			double wt_old_old = getVariable(e).get(k).getScalarVariable(
					VariableName.EQUIVALENT_JUMP, -2);
			double wt_old = getVariable(e).get(k).getScalarVariable(VariableName.EQUIVALENT_JUMP,
					-1);
			MutableDouble wt_alg = getVariable(e).get(k).getScalarVariable(
					VariableName.EQUIVALENT_JUMP_ALGORITHMIC);
			wt_alg.setValue(wt_old + (wt_old - wt_old_old));

			/*
			 * implicit update
			 */
			MutableDouble wt = getVariable(e).get(k)
					.getScalarVariable(VariableName.EQUIVALENT_JUMP);

			// Debugger.watch("wt_old = ", wt_old);

			if (wtilde > wt_old) {
				wt.setValue(wtilde);
			} else {
				wt.setValue(wt_old);
			}

		}

	}

	@Override
	public void FinalizeSolutionStep(Element e) {
		for (Variable v : getVariable(e)) {
			v.save();
		}
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

}
